LANGEVIN MOLECULAR DYNAMICS DERIVED FROM EHRENFEST DYNAMICS 



ANDERS SZEPESSY 

Abstract. Stochastic Langevin molecular dynamics for nuclei is derived from quantum classical molecular 
dynamics, also called Ehrenfest dynamics, at positive temperature, assuming that the initial data for the 
electrons is stochastically perturbed from the ground state, with a large spectral gap, and that the ratio, 
M, of nuclei and electron mass tends to infinity. The Ehrenfest dynamics is approximated by the Langevin 
dynamics with accuracy o{M~^^'^) on bounded time intervals, which makes the 0(M~^i'^) small friction and 
fluctuation terms visible. The initial electron probability distribution is a Gibbs density derived from the 
Liouville equilibrium solution generated by the nuclei acting as a heat bath for the electrons in the Ehrenfest 
Hamiltonian system. The diffusion and friction coefficients in the Langevin equation satisfy the Einstein's 
fluctuation-dissipation relation. The dissipative friction mechanism comes from the evolution of the electron 
ground state, due to slow dynamics of the nuclei. The stochastic perturbation of the ground state can also 
generate a temperature dependent contribution to the drift, depending on the spectral gap of the electron 
eigenvalues. 
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1. Introduction to Ab Initio Molecular Dynamics and the Main Result 

One method to simulate molecular motion is to use quantum classical molecular dynamics (QCMD), also 
called Ehrenfest dynamics, with the Hamiltonian 



(^•^^ ^ ^2^2^ \x,-X,,\ ^ |X„ - X, 

i<k<j<J ' n=lj=l ' ^ "I l<ji<m<JV ' 



= :Hi{X) 



and the energy 

/ ^{xi,...,xj)*H{X)^{xu...,xj)dxi...dx 



XJ 



where the nuclear positions X„ : [0, oo) R"^, n = 1, . . . , and the electron wave function : [0, oo) 
C solve the Ehrenfest dynamics system 

(1.2) MA* = -{^!\dx^H{X')^'), 

=dx„H,{Xt) 

(1.3) idt^^ ^ H{X')^\ 

see [23], [27]. The Hamiltonian is composed of the kinetic energy of the electrons, the electron-electron 
repulsion, the electron- nuclei attraction, and the repulsion of nuclei (with charge Zn), in the Hartree atomic 
units where the electron mass, electron charge, reduced Planck constant, and the Coulomb force constant 
(47reo)~^ all are one. The mass of the nuclei, which are much greater than one (electron mass), are the 
diagonal elements in the diagonal matrix M . The Ehrenfest dynamics (|1.2m.3[) is a Hamiltonian system with 
the Hamiltonian 

M~W + ]^{^,H^>) =:He 

in the real variable ((AT, *'"), (p, **)) {rE,PE), where 2p = MX and ^ =: W + Therefore, the 

Ehrenfest dynamics conserves this energy |pp + 2^^ (4', H^) . The superscript t denotes the time variable 
A* := X{t), in the fast electron dynamics time scale. 

The Ehrenfest dynamics can be derived from the time-independent Schrodinger equation 

1 ^ 

(1-4) (i/-—^ Ax J $ = 



=-Hs 

see [24], and from the time-dependent Schrodinger equation, cf. [2l [23]. The wave function depends the 
electron and nuclei coordinates and on discrete spin states 

(1.5) $(a;\cri, . . . ,a;"',crj, Ai,i;i, . . . , AatjEat), 

which effect the solutions space: each electron spin aj can take two different values and each nucleus can be 
in a finite set of spin states S„; the Pauli exclusion principle restricts the solutions space to wave functions 
satisfying the antisymmetry/symmetry 

. . ,x^ ,aj, . . .,x'',ak, • • •) = ■■,x'',ak,...,x',aj,...) for any 1 < j, fc < J 

and for any pair of nuclei n and to, with A nucleons and the same number of protons and neutrons, 

$(..., X"\ S„, . . . , A", a„, . . .) = (-1)^$(. . . , A", S„, . . . , A™, E™, . . .) 

cf. [3]. We simplify the notation by writing ^{x, A) instead of the more complete (|1.5p . since the Hamiltonian 
H does not depend on the spin of each particle. We use periodicity in the A-domain and set T := M/(LZ) 
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for a cell size L. The Ehrcnfcst solutfon is shown, in |24j . to approximate the full Schrodinger solution with 
error bounded by 0{A'I^^), for the micro canonical ensemble and for equilibrium observables in the canonical 
ensemble, in the case without caustics and with electron eigenvalues having a large spectral gap. The work 
|26| describes with examples some of the weaknesses and strengths of the Ehrenfest approximation. 

The Ehrenfest dynamics can be further coarse-grained by assuming the electron wave function is in its 
ground state - this so called Born-Oppenheimer approximation leads to 0{M~^) accurate approximation of 
observables of the time-independent Schrodinger equation, for the micro canonical ensemble in the case of 
a spectral gap and no caustics, see [21]. The work [21] also shows that Schrodinger equilibrium observables 
in the canonical ensemble, not including correlation in time can be 0{M~^) accurately approximated by 
Langevin dynamics with the Born-Oppenheimer drift and any friction coefficient proportional to the diffusion 
coefficient. 

The purpose here is to study observables including time-correlation in a bounded time interval to precisely 
determine the friction/diffusion coefficient for the Langevin dynamics: Theorem 11.11 shows that Ehrenfest dy- 
namics is approximated by the Langevin dynamics for observables including time-correlation, in the canonical 
ensemble when the initial data for the electrons is a temperature dependent stochastic perturbation of its 
ground state; the accuracy is o(Af ~^/^) on bounded time intervals (in the slow time scale of nuclei motion), 
which makes the ©(A/^^/^) small friction and fluctuation terms visible. The stochastic perturbation of the 
ground state can also generate a temperature dependent correction to the drift, depending on the spectral 
gap of the electron eigenvalues. The next section explains the stochastic setting and the Born-Oppenheimer 
approximation. The main idea in the work is the formulation of a classical Gibbs equilibrium distribu- 
tion of eigenstates, motivated by nuclei acting as heat bath for the electrons in the quantum Ehrenfest and 
Schrodinger Hamiltonian systems. The error analysis uses the residual in the Kolmogorov equation of the 
Langevin stochastics, evaluated along the Ehrenfest dynamics on sufficiently large time steps, to have the fluc- 
tuations behaving as Brownian motion increments; here a critical idea is to use that the Ehrenfest solutions 
are the eigenfunctions of the logarithm of the solution operator. The results are based on the Hamiltonian 
—2"^ 'Y^- /S.xj + Hi with any potential Hj that is smooth in X, e.g. a regularized version of the Coulomb 
potential. 

The particles with coordinates x, in ()1.2|) - ()1.3p . can also be interpreted as a heat bath of lighter particles 
consisting of both nuclei and electrons, i.e. not necessarily only of electrons, so that the Langevin equation 
(|1.15p also describes approximately the dynamics X of heavy so called Brownian particles. This section ends 
with some background on deriving the Langevin equation for heavy particles in a heat bath. Theorem 1 1.1 1 is 
relevant for the central problem in statistical mechanics to show that Hamiltonian dynamics of heavy particles, 
coupled to a heat bath of many lighter particles with random initial data, can be approximately described 
by Langevin's equation, as motivated by the pioneering work [5].|20| and continued with more precise heat 
bath models, based on harmonic interactions, in [T^l [H] [IB]. More recently these models of a heavy, so 
called Brownian, particle coupled to a heat bath are also used for numerical analysis studies related to coarse- 
graining in molecular dynamics and weak convergence analysis |251 119j [15j , for strong convergence analysis 
[1] , and for computational studies on nonlinear heat bath models [3j [18] . Langevin's equation has also been 
derived from a heavy particle colliding with an ideal gas heat bath, where the initial light particle positions 
are modeled by a Poisson point process and initial particle velocities are independent Maxwell distributed; 
the heavy particle collides elastically with the ideal gas particles and moves uniformly in between, see [9l[8]. 

The main inspiration of this paper is [28j, where Zwanzig derives a generalized Langevin equation for a 
(Brownian) particle coupled to a heat bath particle system through a harmonic interaction potential with 
a Debye heat bath frequency distribution and an ad hoc coupling; similar assumptions on deterministic or 
stochastic frequency distributions and coupling are used in the above cited references. This work extends 
the ideas in [28] by using the ab initio Ehrenfest dynamics (|1.2[) for nuclei and (|1.3p for electrons (or (|1.2p 
for heavy nuclei coupled to a heat bath of lighter nuclei and all electrons modeled by a Schrodinger equation 
(|1.3p ). Other differences are that: 
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• the residual of the Langevin Kolmogorov equation along the Ehrenfest dynamics is used to derive an 
error estimate, instead of using an explicit solution, and 

• the slow nuclei dynamics compared to the fast electron dynamics is exploited, in integration by parts 
of highly oscillatory integrands based on the Ehrenfest dynamics solution operator, to find a proper 
Langevin equation, avoiding the integral coupling in the generalized Langevin equation. 

The two central ideas in deriving Langevin dynamics from coupling to a heat bath - to find the friction 
mechanism in the heavy particle coupling to the dynamics of the lighter particles and to find the diffusion 
from fluctuations in initial data of the light particles - were already the basis in [TH |TT] and [5S] . 

1.1. Perturbation of the Born-Oppenheimer Approximation. At low temperature we expect that the 
fast electron dynamics, compared to the slower nuclei in the Ehrenfest dynamics, yields an electron wave 
solution that is almost in its ground state ^fg, which solves the electron eigenvalue problem 

(1.6) H-^a = Ao*o, 

and is normalized (\E'o, ^'o) = 1; here Aq — Ao(^*) is the smallest eigenvalue oi H = H{X*) in L?'{dx). The 
funct 



ion 



4':=exp(-i/ \Q{X'')ds)^o{t,x]X*) 



satisfies 



so that if the nuclei do not move, the wave function ^ solves the time dependent Schrodinger equation; 
and if they move slowly, i.e. the L^-norm H'l'oll is small, the function ^ is an approximate solution to the 
Ehrenfest dynamics. The approximation (X*,\I/), with 5" replacing ^' in (jl.2[) . is called Born-Oppenheimer 
molecular dynamics [23l [14] and it approximates observables of the time-independent Schrodinger equation 
for the electron-nuclei system with accuracy 0(M~^), when there is a spectral gap and no caustics, see [21]. 
The aim of this work is to study molecular dynamics when the electron states ^ are randomly perturbed 
from the ground state, with a Gibbs distribution at positive temperature. To study the perturbation of the 
ground state we first make the transformation 

(1.7) ^{t,x) = e-'^o^°i^'')d'-^(t,x) 
which implies that ^' solves the Schrodinger equation 

(1.8) idt^' = {H{X*) - MX')) ^\ 

^ V ' 

with the translated Hamiltonian i7, and then we make the Ansatz 

= ^'o + * 

and expect : [0, oo) x R'^'^ ^ C to be small. The Ansatz implies that ^ solves 

(1.9) i9t^ = ff(X*)*-i9t*o, 
which yields the solution representation 

t 



^^Stfi^'^^i St.s ds, 



with the solution operator S defined by 

(1.10) St,sip' 

for a solution 

idt(p* = H'ifi^ t > s. 
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The first term in the representation depends only on the initial data and the second term depends only 
on the residual • This splitting, inserted into the equation (|1.2p for the nuclei, eliminates formally the 
electrons and generates fluctuations, from stochastic initial data ^I''^, and friction, through the coupling to 

1.2. Stochastic Electron Initial Data. The next step, to determine the stochastic initial data for ^P, 
requires some additional modeling. Inspired by the study of a classical heat bath of harmonic oscillators in 
|28j , we will sample the initial data for ^ randomly from a probability density given by an equilibrium solution 
f (i.e. dtf = 0) of the Liouville equation dtf + dpj^HEdrEf — c^rE-fffi^pg/ = 0, to the Ehrenfest dynamics 
(|1.2[) and (|1.8p . There are many such equilibrium solutions but there is only one where the momenta pj 
are independent and that is the standard canonical ensemble for Hamiltonian dynamics, i.e. f{rE,PE) = 
cexp (—He/T), where the positive parameter T is the temperature, in units of the Boltzmann constant, 
and c is a normalization constant. There are other equilibrium solutions, e.g. / — h{HE) for differentiable 
functions h. To find a unique equilibrium solution to sample 4' from, we may think of the nuclei as the 
heat bath for the electrons. The equilibrium distribution of the nuclei seems simpler to understand than 
the equilibrium of electrons: in a statistical mechanics view, the probability of finding the heat bath with 
the energy He '■= M~^\p\'^ + {'i/ , H{X)'^) (for given is proportional to the Gibbs-Boltzmann factor 
exp{—HE/T), cf. [To], [17]; an alternative motivation, based on the conclusion of this work (in a somewhat 
circular argument), is that the nuclei can be approximated by classical Langevin dynamics which has its 
unique invariant density exp ( — M~^\p\'^/2 — Ao(X)) /T as the marginal distribution of exp{—HE/T) when 
integrating over all the electron states Since the energy is conserved for the Ehrenfest dynamics - now 
viewed with the electrons as the primary systems coupled to the heat bath of nuclei - the probability of 
finding the electrons in a certain configuration 5* is the same as finding the nuclei in a state with energy 
He, which is proportional to exp(— i/g/T). This conclusion, that the probability to have an electron wave 
function ^ is proportional to exp{—HE/T)d'i'^d'^'^ is our motivation to sample the initial data for 5* from 
the conditioned density generated be exp {— , H'i>) / {2T)): since we seek data for the electrons, we use the 
probability distribution for conditioned on {X,p). 

We compare in Section 15.21 our model of initial data with a more standard model of initial data, having 
given probabilities to be in mixed states. It turns out that Einstein's fluctuation-dissipation result does not 
hold for the canonical distribution, which is a motivation for our focus on the initial data sampled from 
the Liouville equilibrium solution above, where the Einstein fluctuation-dissipation results holds. To sample 
from such an equilibrium density is standard in classical Hamiltonian statistical mechanics but it seems non 
standard for Ehrenfest quantum dynamics. The Gibbs density for the Ehrenfest dynamics approximates 
0{M~^) accurately the Gibbs density for the Hamiltonian dynamics of the time-independent Schrodinger 
equation derived in '24], in the case no caustics are present. 

1.2.1. Slow Nuclear Dynamics. The large mass, M ^ 1, in the dynamics of the nuclei 

= -(**,ajfi/(x*)^'*) 

implies that position increments |X* — = 0{M^^^^) are asymptotically negligible for bounded time 
intervals t = 0{1). It is therefore necessary to follow the nuclei for long time t > M^^^ to see non trivial 
dynamics and then it is convenient to change the time scale by introducing t —: \[Mt f, X'^ and 
:— \lMur to obtain 

fl 111 dr^ - ^ ' 

' ' ^ ^-{^^,dxH{X^)W), 



which means that we want to understand electron wave solutions to (jl.3p up to fast time of order 
t = f ~ A/1/2. 
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To describe the stochastic initial data, we use a spectral decomposition defined in the next section. 

1.2.2. The Initial Data and the Energy Cut-off. Reference shows that observables of the Ehrenfest dy- 
namics approximates observables of the time-independent Schrodinger equation with error 0{M^^). The 
solution of the time-independent Schrodinger equation can therefore be characterized by data for the Ehren- 
fest dynamics. We choose these stochastic initial data sampled from the Gibbs-Boltzmann equilibrium density 
as follows. 

Introduce the normalized electron eigenvectors : j = 0, 1, 2 . . . , J} of the Hamiltonian H{X) : 

iJ^(T^'^) ^ H~^{T^'^) with its eigenvalues {Xj{X) : j = 0,1,2 . . . , J} in increasing order satisfying 
lf{X)ipj{X) — Xj{X)'ijjj{X), for given nuclei position X. We let the initial data for the initial electron 
wave function '5° in the Ehrenfest dynamics (|1.8|) be a stochastic perturbation (|1.9|) of the ground state with 
initial data 

(1.12) = J2^,i.,iX') 

where the complex coefficients 7j = 7j + *7j, for j — 1, ... J , are all independent, with independent normally 

distributed real and imaginary parts of mean zero and variance T/Xj{X^) and 70 := (1 — J2'j=i lljl'^V^^- 
We consider the case when there is a large spectral gap around the ground state energy - more precisely 
we assume that (jl.l4p or (|1.18l) holds, which implies that is a small perturbation of the ground state 
7of/'o(X''). Although the perturbation X]n>o T"^" small its expected energy 

n,m>0 n>0 

is large and we will se that the small perturbation in fact influences the dynamics to order 0{M^^^^). The 
cut-off state J is chosen such that 

(1.13) J^oo 

(slowly) as M ^ 00, to include all eigenstates in the limit. Section [5] presents motivations of the normal 
distribution from the canonical Gibbs-Boltzmann distribution. 

1.3. The Main Result. Let denote the standard Brownian process (at time r) in M"^^ with independent 
components. To simplify the notation we assume that all nuclei have the same mass M 1 — electron mass; 
this can easily be extended to varying nuclear masses much larger than the electron mass. We apply the 
notation ip{x,X) = 0{M~°') also for complex valued functions, meaning that |-0(a;,^)| = 0{M~°') holds 
uniformly in x and X. 

We will use the following .spectral gap condition for the electron eigenvalues 

n 14^ T,i=i = 0(1) and Er=j(V'n, ^x-ff^-n) = o(l) 

We also assume that the electron eigenvalues do not cross, i.e. |A„ — Am| > c for a positive constant c along 
solution paths and n ^ m. The conditions are motivated and used as follows: 

The first condition will imply that the Ehrenfest wave function is a small perturbation of the ground state, 
the temperature is small compared to the large gap Ai and the integral kernel for the friction and diffusion is 
well approximated by a point mass in (|4.2p . Since the eigenvalues may grow slowly for a Schrodinger problem, 
the first condition typically requires J 00 slowly, i.e. the number of exited states J is small compared to 
the gaps: the reason to have J — > 00 is to establish ^'^^ji'^PmdxH^ijn) = o(l) in (|4.3p . which has a high 
convergence rate, since dxH is smooth and ipn becomes highly oscillatory, compatible with J tending to 
infinity slowly. 
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The second condition means that the electron eigenvalues are almost parallel to the ground state eigenvalue, 
as a function of X , which implies that the ground state Born-Oppenheimer drift is a good approximation for 
all states. When this condition is relaxed in (jl.lSp . the drift receives a temperature dependent correction. 
We have dx^n — {''Pm dxH-ipn) which decays fast with growing eigenvalues, since dxH is smooth and ipn 
become highly oscillatory as n oo. 

Let K he a positive parameter; the stochastic Langevin dynamics 

dX^ = p^dT 

dp"- = -dx\o{X^)dT - M-'^/^Kp'dT + y/2TM-'^/'^KdW^ 
has the unique invariant probability density 

^-(pop/2+Xo(X))/Tj:p^X 

cf. [5] . The work |24j shows that Langevin dynamics approximates equilibrium observables of the Schrodinger 
equation and Ehrenfest dynamics with accuracy o(M~^/^) provided (|1.14p holds, in the case where the ob- 
servable does not depend on time and it is enough to determine an integral with respect to the invariant 
measure. Here we study Langevin dynamics approximating Ehrenfest dynamics on a bounded time interval. 
This determines the otherwise undetermined friction and diffusion matrix K, and makes it possible to deter- 
mine observables that depend on the correlation of the positions at different time in the canonical ensemble 
o(M~^/^) accurately. 

We use the notation f = M^/^r for the transformation from slow to fast time, the eigenvalue Aq denotes 
the ground state electron energy (|1.6p of H with normalized ground state \l/o and Xj are the translated 
eigenvalues of iJ — Aq defined in (jl.ip . Define also the rank one dissipation and fluctuation matrix 

K{X) ■.= 2{dx^o{X),dx^o{X)), 
K^'\X) - V2|ax*or'(5x*oW,9x*o(^)>. 

Theorem 1.1. Assume that the spectral gap condition \1.1A\ holds, the temperature is bounded T = 0(1), 
and the electron initial data is given by (jl.l2p . then the ltd Langevin dynamics 

, , dX^ =-dx\o{X^)dT-M-^/^K{X^)X''dT + V2jW^K^^^{X^)dW^, 

^ ^ o<r<r, 

approximates Ehrenfest dynamics with accuracy 

(1.16) \E[g{X^,p^)^g{X^,X^)]\^o{M-^/^), 

for any bounded function g : MP^ x M.^^ — > R, provided the Langevin value function 

u{X, p,a) := E[g{X'^, X^) \ X" = X, X" = p] 
has bounded derivatives of order one to four 

(1.17) / s\rp \Du{X,p,a)\iida ^0{l), D ^ dr, dx , dp, . . . , dpppp. 

Jo {X,p) 

The approximation result uses a non interacting particle, with given velocity equal to one and position 
coordinate Xq — r, that acts as the time coordinate, so that e.g. transport coefficients as diffusion can be 
studied. 

The small dissipation term M~^^^Kp and small diffusion term are visible in the convergence rate o{M~-^/^). 
The approximation of Schrodinger observables in equilibrium, in |24| . is based on long time and determines 
also the diffusion matrix to be equal to the dissipation matrix. That is, the Langevin dynamics (|1.15|) satisfies 
the Einstein's fluctuation-dissipation result: the square of the diffusion coefficient is the dissipation coefficient 
times twice the temperature. 
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The combination of the cquiUbrium resuh from |24j . that Langevin dynamics approximates the equihb- 
rium observables of the Ehrenfest and Schrodinger equation, and the approximation of Ehrenfest correlation 
observables over bounded time intervals 0(1) here shows that this Langevin dynamics is an accurate approx- 
imation of Ehrenfest dynamics for both short and long time. We cannot from this conclude that Langevin 
dynamics is also a good approximation of the Schrodinger equation for short time correlations but conjecture 
this is true, considering the qualitative similarity of the Ehrenfest dynamics and the Schrodinger dynamics 
as formulated in 24J, where H is essentially replaced hy H — ^j^(2M)~^ Aj^^ . 

Theorem 1.2. If we instead of the gap condition (|1.14p assume the weaker condition 

n 18) ^Li ~^n' = 0(1) and EZA^n^dxHi^n) = o(M-V2) 

E„>o(|5^.A„|A-i)2 =o{M-^/^)andj:n>o\dx,Xn\X-'^0{l), 

then the result in Theorem \l.l \ holds, provided the drift dx\) is replaced hy 

dx{MX) + |trace^ \ogH[X)), 

where trace^ is the trace in the orthogonal complement of 4*0 . 

2. The Force on Nuclei 

It is convenient to split the force on nuclei in (|1.2p into two parts, using the definition of in (|1.7p . 

(*,9xi?(^)^) =(*, 9x^(^)4') 

- {^,dxMX)^) + {^,dx{H{X) - Xo{X))^) 
^dxXo{X) + {^,dxH{X)^). 

With the Ansatz = ^'o + the second term in the nuclear force becomes 

(*,ax-ff(X)*) = (*o,ax^(^)*o) 

+ {^,dxH{X}^o) + {^o,dxH{X)^) 
+ {^,dxH{X)^). 

Use H'^o — and consequently 

(2.1) dxH^o + Hdx'fo^O 

to obtain for the first term 

{^o,dxHiX)-^o) =-{^o,HiX)dx'^o) 
= -(£*o,5x^'o) =0. 

=0 

Let 3? denote the real part. The second terms are 

2^ dxH{X)^o) + (*o, dxHiX)-^) = 2^{^>, dxH{X)-^a) 

^ ■ ' ^~2^{^,H{X)dx^o) 



which, by the dynamics (|1.9p , generates a fluctuation term 

(2.3) -23fi(5t,o*(0),iJ(X*)9x*^) 
and a friction term 

(2.4) -23fj/ {St,sdx%X',H{X')dx'<i^l)ds. 

Jo 

In Section l4?2l we show that the third term, {'^,dxH'^), is negligible small for spectral gap condition (|1.14p 
and yields the correction r9j5s:trace^if/2 for the weaker gap condition (|1.18p . The asymptotic estimates of 
forces are explained more precisely in Sections 14.11 and 14.21 
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Remark 2.1 (Quantum initial data). Assume, instead of (|1.12p . that is a pure electron eigenstate e'^^^tpj 
with Gibbs-Bohzmann probabihty 

(2.5) ,^.:=e-^^/^(^e-^^/^)-i 

e 

(and independent random phase shifts aj uniformly distributed on [0, 27r]) for j — 0, . . . , J, and write 
>o ^i^i which has the covariance E[(7j)*7fc] — q^S^k- Let ^'o(-''^) — 7o V'o(-'^)- Then the fluctuations 
are very different from the case in Theorem 11.11 since E[{^{^° , Hdx^o))'^] is zero due to ^[7*7^ 7970] = 
for j, k > 0. We assume in (|1.12l) that the electron-nuclei system is in a pure eigenstate of the full Schrodinger 
equation p.4p and that does not mean that the electrons have to be in eigenstates of the electron operator 
H for fixed nuclei positions. 



3. Approximation Error of Observables 
3.1. The Dynamics. The Ehrenfest dynamics can be written in the slow time scale 

p- ^-dxXo-{r,dxH{X^W) 



with the Hamiltonian 



using the coordinates {X,^r',P,^i) and -0 := -0^ + iipi — {2/My^^(j}. We shall approximate the Ehrenfest 
dynamics by (X^.p^) defined from the Ito-Langevin dynamics 

^ ' PL - -dxXoiXL) ~ M-^/^K[Xl)pl + {2TM-^/^f/^K^/\XL)W, 

for gap condition (|1.14p . and with dx{^o + trace^i/) replacing the drift 9xAo for gap condition (jl.lSp . To 
simplify the analysis of the coupling between {X,p) and the thermal fluctuations, introduce the electron wave 
functions {V'n} resembling individual states 

(3.2) -^ijl^HiX^Wn. 0°=0„(X°), 

which implies 

J 

n=0 

where 7„ are independent normal distributed with mean zero and variance T/\n{X^). The Schrodinger 
dynamics (|3.2p shows that {0,^ | n = 0, 1, . . .} forms an orthogonal set 

,3 3^ |:(0„, 0m) = (-lMl/27J0„, 0„) + (0„, -iM^/^HlPm) 

^ ' ' = (0„, iAfi/2iJ0;„) + (^„, -zMV2iJv:„) ^ 0, 

since the initial data {^p^ | 71 = 0, 1, . . .} is orthogonal. 
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3.2. An Error Representation. Define for the given observable g : M'^^ x M'^^ M and the Langevin 
dynamics in (|3.1|) . the value function 

u{y,T) ■.= E[g{Xj,pl) \ {Xl,pl)=y], 

which solves the Kolmogorov backward equation 

dru+podxu - {dxM^) + M-^''^K{X)p) o dpU 
(3.4) +Y,^M-^'^TK{X)dp^p^u^Q r<T 

where t; o u; is the usual Euclidean scalar product in R^^. The goal is to analyze the error ]E[g(X^,p^) — 
of the heavy nuclei particles in the Ehrenfest dynamics approximated by the Langevin dynamics. 
This error can be written as the residual of the Langevin Kolmogorov solution (|3.4|) along the Ehrenfest 
dynamics {X'^ ,p'^) 

ng^X^^p-") - g{Xl,pl)] = nuiX^^p", T) - u{XlJl 0)] 

= E[u(X^,p^,T)-u(X0,pO,0)] 

-nu{Xl,ptQ)-u{X^yM- 
. ' 

=0 

We assume that initial data for the Langevin dynamics is the same as for the Ehrenfest dynamics. Therefore 
the last term is zero. The first term in the right hand side will be written as an integral over time and the 
assumption p.l7p makes the integral bounded. 

We will see that the fluctuations in Ehrenfest dynamics behaves approximately as the Wiener process 
fluctuations in the Langevin dynamics only on sufficiently large time intervals At. Therefore we per- 
form the analysis on such large time steps and introduce for time steps At the notation X(jAt) := X\ 
dX'^ /dT\r=jAT = p(j At) —: pP and Tj := jAt. Telescoping cancelation implies 



w(XO,pO,0)] = E,E 



u(X^+i,p^+i,T,+i)-u(X^p^T,) 

u(XJ + l,p^ + l,T, + l) - u(X^p^Tj + l) 



+ u(X^p^T, + i)-u(X^p^T,) 



and Taylor expansion, with the position and velocity increments 

AX := X^+i - X^ = 12^^ p'^da, 



shows that 



Ap :=pJ"+i -p^ = -J^^+'{r,dxHiX'')r)da, 
E[u(XT,j5T,T)-u(X0,pO,0)] 

= e[9.u(X^pJ,t,+i)At + dMX\p',T,+i)AX 



+dpu{X^ ,p> ,T,+^)Ap + ^dppu{X^,p^, T,+i)(Ap)^ 



+ error termsi, 



for some point (X^ jp') on the line between (X^ ,p^) and here the error terms are 

error termsi ■■=Ej^ [drMX^ , p' , fj ) ( At)2/2 + d^^X^lP , r^+i ) (AX)2/2 
+d^pu{X\p^,Tj+i)AXAp 
for some t, between Tj and t^+i, using the notation 



wdyU 



d^yuAXAY = ^^.fc d^^y.uAXjAYk, 



dyyUl 



j '-'VjVj ' 
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Apply the Kolmogorov equation to eliminate the drU term 

in the error representation to obtain 

E[u(X^,p^,T)-u(X0,pO,0)] 
= Y.j^[dxu{X^ ,Tj+i){^X - pi ^t) 

(3.5) +9p?/(X^p^r,+i)(Ap - (-9xAo(XJ)-M-i/2x(XJy)Ar) 

+i(appu(x^p^T,+i)(Ap)2 - dppu{xi ,p> ,T,+i)j^k(xi)AT) 

+crror tcrmsi. 

These expected values are determined in Sections I4.1H4.3I Section 14.41 estimates the three error terms 

(3.6) ^[drru{l^T)% ¥.[d^Ml^X)\ E[9^„uAXAp]. 
We have 

Api := - /; {dxMXD + M-^I^K{Xi)pl) ds + /J {M-^I^TK{Xi)Y'^ 

and 

Ap = - /; (9xAo(X«) + E„,,„ inl*m{rn,dxH^rj) ds 

(3.7) = -/; (axAo(X^) + |7oP(V;o,axi/^o) + 2E„^o3?(7„^„,9xi?V;o)) ds 

In the next section we show that the last integral with the quadratic terms is negligible small, due to our 
assumption X]n>o^n^ ^ o(Af^^/^) in condition (|1.14p . and yields the correction of the drift for condition 
(fTTSl). 



Introduce the orthogonal decomposition 

(3.8) ^„ = ^„©^^, 

based on the electron eigenvectors i/jn satisfying H-ipn = Ki'4'n- The Born-Oppenheimer approximation (|4.15p 
shows that 

(3.9) ^i[t)^0{M'''^) 

and since H^jq = we have by (|4.14p the following representation in the slow time scale 

i,^{T) - ^.,o^(}(0) - i I ~SrA^o)\a)da, 



where -ij^Q = ('^o, 5'o)^'offi'0o the orthogonal decomposition and S is the solution operator Sr,a'4'Q{o') — iPq{t) 
in the slow time scale. The ground state relation Htpo — implies dxH^o + Hdxi^o — 0, which combined 
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with = 0(M-i/2) yield the main friction term 




=0 =o(M-i/2) 

(3.10) = -2di{'^^,Hdx^o) + o{M-^'^) 




+o(Af-i/2). 

4. Estimates of the Friction, Fluctuation and Coupling 

In the error estimate there are two main terms: 

, , a time discretization error 0(Ar), and 

^ ' ' a couphng error ~X^^O{{MAt)-^) + o{M-^/^). 

These two errors combine to the bound 0(\i ^^^Af^^/^)+o(A/^^/^), for the optimal choice At = ^^^M^^/^, 
to obtain the error estimate o(Af ^^/^) in p.l6|) . 

The error representation ()3.5|) consists of a X-drift term with the factor dxu, a p-drift term with the factor 
dpU, a diffusion term with the factor dppU and the three error terms (13. 6p analyzed in the following sections. 

4.1. The Priction Term. This section evaluates the friction term (I3.10p using the solution operator ()1.10|) . 
its time derivative and integration by parts (in the slow time scale) to remove the factor H: study the first 
term in the right hand side of (|3.10p 

-25R(/q iStM ds, H*dx%) 

= ~2M-y^^{j;^f^{SUHT'n)ds,H'dx^h) 

+2M-^/^^{J^St,s4-AiHT^n)d-^^H'9x%) 

- -2M-^^^^{St,t{H')-^%,H'dx'i'h) 
(4.2) +2M-V2gfi(^,_o(7fO)-i^o^ H'dx%) 

+2M-^/^^{J^St.s£{{Hn-'n)ds,H'd^q.t^) 

^-2M-'/^j:^{dx,Kdxn)X, 

+2M''^/^^{St,o{H°)-'^^°o^ H*dx%) 

+2M-y^^{J^St,ff^{{HT^n)ds,H'dx%). 

Here we use that the null space of H{X) is spanned by 'i'Q{X) which is orthogonal to ^o(-''^) = po dx^o{X), 
so that H~^{X)^o{X) = X]n>o '^""^^'"(V'nj ^o) is bounded by the assumption of a spectral gap (I1.14p or 
(|1.18p . The first term in the right hand side yields the main friction term 2Af~^/^ Jr^^^ {dx'^o, dx'^njpds. 

The integral over one time step of the second term above can also be integrated by parts as above to gain 
a power of Af loosing one power of (Ar)~^ and gaining one power of Aj^^ = oil) from i?^^, contributing 
in total to the coupling error 0(A]"^ Af~^ At~^). 

The last term in the right hand side above can be integrated by parts again to give a negligible error 
contribution of order 0{M~^). 

4.2. The Fluctuation-Dissipation Property. This section derives the fluctuation dissipation property 
based on the representation (p ~ with independent normal distributed random variables 7„ and the 
orthonormal set {V'™ | n = 1, . . .}, orthogonal to V'o- We use the notation a ~ &, meaning a = & + o(l), at the 
following five steps in the calculation below: 
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1. in the first step the error term comes from being slightly dependent on 7„ - this coupling 
yields a small error term estimated in Section 14.31 

2. the second step uses the spectral gap condition (I1.14p to relate Aj'j and A^j, explained more below, 

3. the third step uses that ^„ = + C'(Af"^/^) in (j3.8p and ()3.9|) to replace a factor of with ^ 

4. the fourth step applies assumption (|1.13p . that J oo, and the orthogonal set {V'n}5^Li in (|3.3p to 
have a sum of eigenstates forming a basis in the L^(T'^'') orthogonal complement of tpo, and 

5. the fifth step uses that 'tpo is orthogonal to dx'^o and -00 — "00; 



+25iEtiEi 



__,{Sr,.i^Z,Hdx%){^?^,Hdxn 




(4.3) 



Integration J^" 



^mE'U{Sr,.4'Z,Hdxn){dxn.H4,^)^ 

^ 4.mj:i^,{Sr,.i;!^,Hdx%){dxn,4'Z) 

^ 4T3?Er=i (C, {Sr,.rHdxn) {dxn. C> 

^^Tnj:n=oirnASr.arHdx%){dx^^0.rn) 

^4:m{Sr,adxn.H^dxn)- 
. da gives 2r times the friction term we have in the left hand side of 



im{Sr,,^,dx'f^, H^-dx%-)do 



.-At 



without the factor p — X, with the difference that the integration is over a time interval of length At instead 
oft. Since this is also integrated over a time interval with respect to r, we conclude that the main contribution 
and the error terms are the same as for the friction case. 

We see that this fluctuation-dissipation result requires precisely that the variance of 7„ is r/A„ (to leading 
order), that E[7„7m] is zero for n ^ m and that {ipn}n>o forms an orthonormal basis. In the second step we 
used that 



(4.4) 

where the error term 
(4.5) 



1 



1 



-(1 



{X° - X-) O dxXn 



A„(XO) A„(X" 

^ (X" - X-) O dx~Xn 

n>0 



At] 



)• 



0(1) 



is negligible by the spectral gap assumption (|1.14[) or (jl.lE 
4.3. Couplings. The p— drift has a fluctuation part 



(4.6) 



Ap 



This integral would have vanishing expected value if X'^ and tpn would be independent of 7„, since E[7„] — 0. 
In this section we show that the dependence between {X,p,(j)) and 7„ is asymptotically negligible, so that 



14 



ANDERS SZEPESSY 



the expectation of this fluctuation part of the drift becomes small 

(4.7) E[apM(X^",/",f„+i)Ap] = Ato(A/-i/2). 
The estimate is subtle in the sense that we expect that Aj5 behaves as 

and we will need to use that dpU is evaluated in the point p"^" and not p'^"^'^'^, since e.g. a factor p in dpu{X, p, •) 
would give the coupling error E[Ap(g) Aj5] = 2TM-^/'^K^/'^At by (gSl), which is too large compared to ([i??)) . 
To establish the improved estimate (|4.7p we use the oscillatory behavior of ip^ in integration by parts, as in 
the two error terms of (j4.2p . which distinguishes between the arguments p'^"+'^'^ and p"^" (since r will always 
be after t„). In this sense the product of the fluctuation term (|4.6p and 9pM in (|3.5|) generates one of the 
most important terms - the one corresponding to ltd martingales. 

Define the first variation d-y^{X,p,(j)) =: {X'j^,p'j^,(f)'^)^ which satisfies the linearized Ehrenfest system 

X' ^p', 

p' = -X' o dxxMX) - 2^{i^.,dxH{X)4>) ~ (0, X' o dxxH<t>) 

(4.8) _ -25R(E„7nVi;,5xi?0), 

j^i,'^ =Hi^'^ + X'odxHi''n, 
(X,p,</))'(0) -0, 

where (f) = ^„ ^nipn- Duhamel's principle shows that 

iX.P,mT) =-2/;/(r,a)3?(Vi„,Sx^(X)<^)(a)da 
^-2/;/(r,a)(C,ff5x*^)da 

and / is the linear solution operator for (|4.8p . which is bounded 0(1). Use this and the first variation to 
write the coupling as 

(4.10) E[V(^"",P"",^«+i)Ap]= / E^" / 5d7«dT, 



ri>0 



where the function g is the first variation of dpu{'ipm Hdx'^o) a-nd by the chain rule it is a function of {X,p, (/))- 
derivatives of dpu{ipm Hdx'^a) and of /. The coupling is then similar to the fluctuation term determined 
in Section [42l with a product (7„V'^, 7J'^9x4'o)(7nV'ni -^'^^^^o) the couphng with dpU; the difference 
that makes this coupling smaller is that r and a are not in the same time interval, so the main term in 
(|4.2p will not appear, instead there will be two terms like the second term in the right hand side of (|4.2p 
containing the solution operator, generating an extra factor o(l) of decay in the case At ^ M^^/^. The 
coupling (7n'0n? 9x{H'^dx'^o)) (7nV'n' H'^dx'^o) introduced from p itself can be handled by noting that 

X\T)=X\Tn)+ p'{s)ds 

where the X'(t„) term behaves as above; the second term now has coupling in the same time interval, but 
here there is an additional factor of At coming from the time integral yielding the total coupling bound 
e)(M-i/2(AT)2) satisfying (jiT]) . 

We can similarly determine the coupling of quadratic terms. The diffusion term 

E[5ppu(X",p",T„+i)Ap® Ap] 

has the main contribution from 

E[appu(X",p",T„+i)Ap® Ap]. 
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To analyze the diffusion term similarly as the drift term, we differentiate (|4.9p with respect to 7m to obtain 
a linear system for the second variation d'y^'y^{X,p^ </>) 



4>)'Lm{r) = -2 / /(r, amil;n.dxH{X)^„,){a)da 
Jo 

and we obtain as in (|4.10|) 

The second variation shows also that the remaining coupling terms 

E[apMAj5] 

E[9ppM(X",p",T„+i)Ap(8 Ap] 
E[appw(X",p", r„+i)Ap Ap] 

E[appw(X",p",r„+i)(Ap® Ap-2ri^i/2M-i/2AT)] 

with 

lnlm{'4'n,dxH^pm)'j ds, 

are bounded by Ato{M~^^^). 

The second term inp uses E[7„7m] = and coupling as follows. The solutions tpn are the eigenfunctions q„ 
of the essentially self-adjoint operator Q, defined by the backward solution operator 5*0,* — e^'*^ ' '^o.t 
V'n(O) = 'S'o,t'0ri(i)) with real eigenvalues A^, explained as follows. Let /(O) — e'^^ ^°'*qn{t). The continuity 
at zero gives /(O) = qn{0) and the limit t 0+ in 

t^^Qn,tqn{t) = K{t)qn{t) 

combined with t^^Qa.t i?" shows that (7„(0) are the eigenfunctions of H'^. We have obtained 

which shows that qn{t) is the solution of the Ehrenfest dynamics along X with initial data ipn{0) and conse- 
quently 

Differentiation of the relation Qo^tqnit) — ^n{t)qn{t) and Q ~ H shows that A„(0) — A„(0). Therefore 



is oscillatory as a function of s and integration by parts in the coupling term 
(4.11) 9^\ln\^\l,n\''{MdxHMV{MdxHMyda 



yields by the gap condition (|1.14p or (jl.lSp 

E[^ /_] = E[ 5] |7„P ^ |7™P] = o(M-i/^), 

n,m n m 

provided the phases A„ — A„i have no critical point, which holds if the eigenvalues A„ do not cross as assumed. 
The first term in p for case (|1.14p can be handled without coupling 



nEn>o\lnnMdxHM\] =nEn>0\ln\'{\{MdxHM\+O{M-'/^))] 

+c 
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The first term in p for case ()1.18|) . The first term gives the leading order contribution 
The gap condition (I1.18|) and (|4.4|) show that 

which estabhshes the leading order correction to the drift of the Langevin equation 

J^E = IrdxY^^ogUx') = iTaxtrace^logff(X*) 

where trace"*" is the trace in the orthogonal complement of and the factor 1/2 comes from 9ppu(Ap)^/2 in 
The coupling for this trace drift has, as in (|4.1ip . the factor 

\ln\H'^~M,dxH^nY{in,dxH^nY 

n>0 

with its expected value bounded by En>o(^''f ^"^n ~ o(M^^/^), which makes the coupling negligible. 

4.4. The Three Error Terms. Our main assumption, that the nuclear force is bounded {'if,dxH{X)'$) = 
0(1), implies directly that the three error terms in the error representation (j3.5p have the bound 

|E[9,,u(At)2]| + \E[dxxuiAXf]\ + \E[dxpuAXAp]\ = ©(Ar^). 

When the mean-field drift for the Langevin dynamics is small, the energy e'^ :— |p£p/2 + \o{X'[) is low; 
for instance at the minimal energy the mean-field drift is zero, which defines the equilibrium points. This 
energy evolves by Ito's formula as 

de^ =pIo dpi + dxXoiXl) o dXl + ZTM-^/^K^NdT 

= M-^/'^K^iTN - \pl\'^/2)dT + (2rM-i/2ir^)i/2p£ o dW 

so that its expected value satisfies 

(4.12) dE[e] ^ ¥.[M-^''^K^CiTN - \pl\'^ /2)\dT. 

To have initial data near equilibrium means that the energy is low (i.e. the mean-field drift is small) and 
(|4.12p shows that the mean energy grows until the kinetic energy |pip/2 reaches the value 3TN. This growth 
takes long time of the order Af-*^/^ in the slow time scale. Therefore the Langevin dynamics stays long time 
near equilibrium, where the mean-field drift is small, if it starts from an equilibrium state where the energy 
is low. 

4.5. The Born-Oppenheimer approximation. The purpose of this section is to study the evolution p.2p 
of 4'n and explain ip:^ — 0{M^^/^) in p.9p for the decomposition V'n = "^n ffi "^n j where is an eigenvector 
of H{X'^), with norm of order one, satisfying H'^ipn — Kii'n for an eigenvalue S M. A more careful study 
of the Born-Oppenheimer approximation is in Let := e'^'^ foKds^t^_ xhis A nsatz is motivated by 
the zero residual 

(4.13) #n - M^/^{H - K)4>n = 
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and the small residual for the eigenvector 

M^'^{H -\n)i>n =0, 

where w = | V'n I ~^ (V'n ,w)ipn®uj'^ denotes the orthogonal decomposition in the eigenfunction direction and its 
orthogonal complement, as in Section We have 

and by 



(4.14) ^;^(t) = ^t,o^A^-^*^M(£fn(f))^'^s, 

where S is the solution operator >S't,o'0n — i'n- Therefore the integration by parts (|4.2p can be applied and 
we obtain 

(4.15) V^(0=O(Af'^/^), 
provided Am — An does not vanish, which we assume. 

5. The Constrained Stochastic Initial Data 

As in models of heat baths we assume that the initial data of the light particles (here the electrons) 
are stochastic, sampled from an equilibrium distribution, see [121 IHl El]- We use the Gibbs-Boltzmann 
distribution proportional to exp (— iJ^/T), which is an equilibrium solution of the Liouville equation of the 
Ehrenfest dynamics, as explained in Section [1.21 Let us now determine precise properties of this distribution 
generated by the Hamiltonian He- Diagonalize the electron operator H{X*) by the normalized eigenvectors 
and eigenvalues {^pj,Xj} 

(vI/,ff(X*)v|/)=^A,|7,f 

where 

(5.1) HiX')^P, ^~X,^j, 

Ao -0, 

with real and imaginary parts 7j =: 7J + 17*. We show in Section l5.ll that the equilibrium density can be 
approximated by the decoupled density 

(5 2) D n,>oe-^^l^^l^/^rf7;rf7jd^rfp 



which implies that {7J , 7] , j > 0} are independent and each 7J and 7* is normally distributed 7™ ~ N{Q,T/Xj) 
for j > 0, m = r, i. Note that K[-fj] = for j > and let '^q be the mean field approximation 

*o(^) ■.^Y.^h,]MX) = E[^o]MX) E[(l - ^ |7,f )i/2] ^o(x), 
j>a j>o 

which requires the upper bound on the temperature, TX]j>o 

Note that our initial state is a pure eigenstate of the full Schrodinger operator with energy E and since 
there are many such states we use the model where the initial data is in any state ^P, for ||^|| = 1, with 
probability weight 

e~'''''"''^/^d^'-d^\ 
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The orthogonal transformation ^' = "^jfjipj shows that this probabihty is approximated by D in (15. 2|) . 
using that the determinant of the matrix of eigenvectors is one, where the approximation error comes from 
neglecting the constraint ||^'|| — 1. The advantage with the initial distribution derived by D is that it generates 
a Langevin equation (which also satisfies the classical fluctuation-dissipation result), while initial data given 
by (|2.5p . or a mixed-state variant of that, lead to different fluctuations, see Remark l2.1l and Remark l5.2l Note 
that if we consider densities that are functions of the Hamiltonian He, it is only the exponential function 
(i.e. the Gibbs distribution) that yields independent electron states in this sense, since independence means 
that the joint distribution must be a product of individual distributions. 



5.1. The Constraint ||^'|| — 1. The distribution (|5.2p is determined by neglecting the constraint ||^'|| = 1. 
Let us now motivate it also with the constraint, asymptotically as M tends to infinity. The constraint can, 
by homogeneity of the norm, be explicitly removed if one instead considers the density generated by 



(5.3) 



■E,>o^jl%lV(TE,>ol7,l') 



for all values of jj. This gives a uniform distribution of 70 (conditioned on 7^, j > 0) on the domain of 70, 
which is arbitrary large. Therefore the coupling sum X]j>o ^he denominator is dominated by the term 

|7op. Normahzing jj, j > by 7^ =: 7^-70 reduces (|5.3p to 

e-2:,>o^.l7,lV(T+TE,>ol7,l')^ 

If we let 7j be independent normal distributed with zero mean and variance aT/Xj, for a certain constant a, 
to be determined, the couphng term 1 -I- X]j>o ITj ^i^) close to the expected value 1 + aT J2j>o ^7^' ^^^'-^ 
its variance X]j>o 3r^Aj^ is o(l), as AI tends to infinity. Therefore the condition 1 + ctT^j^^ Xj^ — a, with 
the positive solution a — 1/(1 — TX]j>o J^^^ requiring the condition 

(5.4) Tj2>^j'<h 

j>0 

implies 

We have J2j>o l^j l^ ~ l7opQ; and the normalization leads to 

E,>ol7,f/(«l7on ^tj:^^,xj\ 

which is the same as choosing 7^ from (|5.2p for j > and let I70P = 1 — X]j>o l7jl^- 
The normalization J2i>o l7il^ = ^ couples the variable 70 with the others 



70 = ji -y^i7ii 



j>0 



At zero temperature jj = for j > and 70 = 1. This normalization only makes sense for 1 — X]j>o l7jl^ — 
which will be violated almost surely with normally distributed variables 7^ . Therefore we make a tiny cut-off 
in the probability density to have 1 — J2j>o l7jl^ *-* ''^ith probability one. In fact the analysis uses only the 
mean, the variance and independence of jj, j > 0. 
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5.2. Entropy and the Standard Canonical Density Distribution. Let qj denote the density of state j in 
the initial data '5". hi the usual setting of a canonical Gibbs-Boltzmann distribution qj = e^'^^^^/^^ e"^^/"^, 
as in (|2.5p , which follows from maximizing the von Neumann entropy defined by — qj log qj , with the 
probability and energy constraints qj = 1 and ^jlj — constant, see |13j . The stochastic model for the 
variable |7jP, measuring in (|5.2p and (|5.ip the probability to be in electron state j, is different: the chi-square 
distribution of Aj|7jp/r contains both the weight to be in electron state j and the spacial density of this 
state, see pi] . 
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